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ABSTRACT 

This paper describes a general investigation of stationary oscillations of galaxies. It 
begins with a linear analysis of modes of oscillation with continuous spectra of real 
frequencies. Such modes are gravitational analogues of the van Kampen modes of 
oscillation in plasmas. The characteristic value problem governing modes of the van 
Kampen type in a galaxy is solved with the aid of a modified version of the matrix 
method of Kalnajs in which the perturbation of the distribution function is expressed 
in terms of generalized functions. In general, there is no characteristic equation gov- 
erning the frequencies in the continuous spectrum. However, isolated frequencies in 
the continuous spectrum do satisfy a characteristic equation which, for stellar systems, 
is a counterpart of the dispersion relation proposed by Vlasov for plasma oscillations. 
The linear analysis also provides a characteristic equation for modes with a discrete 
spectrum of real and/or complex frequencies. The second part of the paper describes 
a perturbation theory for a stationary oscillation of a galaxy with a small but finite 
amplitude. Integrals of the stellar motion are constructed with the aid of canonical 
perturbation theory and used in conjunction with the theorem of Jeans in order to 
specify the density of stars in the six-dimensional phase space. These oscillations are 
slightly nonlinear counterparts of the modes of the van Kampen type, and they are 
stellar-dynamical counterparts of the nonlinear plasma waves described by Bernstein, 
Greene, and Kruskal. Fully nonlinear models of stationary oscillations of galaxies can 
be constructed with the aid of Schwarzschild's numerical method for the solution of 
the fundamental integral equation describing the self-consistency of a stellar system. 

Key words: stellar dynamics - galaxies: kinematics and dynamics: structure - insta- 
bilities - plasmas. 

1 INTRODUCTION 

Stationary oscillations of galaxies are oscillations with constant amplitudes. Relatively little attention has been given to 
stationary modes of oscillation in conventional investigations of small perturbations of galaxies. A major reason for this 
apparent neglect is that the principal investigations of the subject have been based on methods which have their origin 
in Landau's (1946) initial-value formulation of the theory of plasma oscillations. A particular feature of such methods is 
the adoption of the so-called 'Landau contour' or 'Landau prescription' for the integration of singular quantities over the 
momentum space (Fridman & Polyachenko 1984; Palmer 1994). Early investigations of the oscillations and the stability of 
infinite, homogeneous steller systems (Lynden-Bell 1962; Sweet 1963) are based on Landau's formulation of the initial-value 
problem. Landau's prescription also underlies the matrix method formulated by Kalnajs (1977) and used by many others (e.g., 
Polyachenko & Shukhman 1981; Fridman & Polyachenko 1984; Palmer & Papaloizou 1987; Weinberg 1989,1991a,b; Saha 1991; 
Bertin et al. 1994; Palmer 1994; Robijn 1995) for investigations of instabilities in galaxies. Such methods are suitable for the 
investigation of instabilities. However, in place of such initial-value formulations, one needs a normal-mode analysis in order 
to study stationary oscillations of a galaxy within the framework of linear perturbation theory. 

The understanding and expectations that one has regarding small perturbations in stellar systems rest heavily on what 
is known about plasma oscillations. In particular, it is well known (see, e. g., Jackson 1960; Clemmow & Dougherty 1969; 
Ecker 1972; Stix 1992) that the initial value formulation of Landau (1946) and the normal mode analysis of van Kampen 
(1955, 1957) provide mutually complementary descriptions of electrostatic waves in a homogeneous plasma. The van Kampen 
modes with continuous spectra of real frequencies are stationary modes of oscillation. However, those modes are singular 
in the sense that perturbations of the distribution of particles in the six-dimensional phase space of a single particle are 
expressed in terms of generalized functions. It is frequently suggested that such distributions are artificial and would be 
difficult to establish as initial conditions in a plasma. On the other hand, the van Kampen modes are complete (Case 1959), 
and they may be superposed in order to represent more physical initial conditions. As results of dispersion and phase mixing, 
such superpositions evolve consistently with Landau's solution of the initial value problem and exhibit forms of apparently 
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irreversible behavior such as Landau damping. Conversely, when a van Kampen mode is adopted as an initial condition, the 
solution of the initial value problem in accordance with Landau's prescription is that van Kampen mode. Thus, the study of 
the van Kampen modes, the linear modes of stationary oscillation in a plasma, contributes significantly to our understanding 
of plasma oscillations, notwithstanding that such modes might be difficult to excite individually and superpostions of such 
modes would suffer Landau damping. 

The literature of plasma physics shows that van Kampen modes are significant and important in other respects as well. 
The van Kampen modes are the small-amplitude limits of nonlinear plasma waves ('BGK waves') of the kind described by 
Bernstein, Greene, and Kruskal (1957; Bohm and Gross 1949; Jackson 1960; Stix 1992). The nonlinearity of BGK waves 
precludes a principle of superposition. Thus, BGK waves are stationary oscillations of a homogeneous plasma, and they do 
not suffer Landau damping. As the linear limit of a BGK wave, each van Kampen mode must be considered individually. 
The interpretation in terms a BGK wave shows how a van Kampen mode owes its self-consistency to the way in which 
orbits trapped in the electrostatic potential of the wave are populated with particles. The connection also establishes the van 
Kampen modes as a starting point in linear perturbation theory for the study of a certain class of nonlinear waves. 

For stellar systems, this paper describes an investigation of stationary oscillations which are the gravitational counterparts 
of van Kampen modes and BGK waves. Section 2 describes the construction of modes of oscillation and instability of a stellar 
system within the framework of linear perturbation theory and with the aid of a modified version of the matrix method 
of Kalnajs (1977; Fridman & Polyachenko 1984; Palmer 1994). In particular, this formulation describes singular modes 
with continuous spectra of real frequencies. The treatment unites the methods of Kalnajs and van Kampen and generalizes 
both. Section 3 describes the construction of slightly nonlinear oscillations of stellar systems. In order to solve the governing 
equations, we make use of canonical perturbation theory (Goldstein 1980) in order to construct suitable integrals of the motion 
for the stellar orbits, and we apply the theorem of Jeans in order to specify the distribution of stars in the six-dimensional 
phase space. The oscillations derived in this way are slightly nonlinear counterparts of the singular modes of the van Kampen 
type that are described in Section 2. As in the case of BGK waves, such stationary oscillations of a galaxy owe their self- 
consistency to the way in which resonant orbits are populated with stars. The nonlinear treatment imposes certain constraints 
on the frequencies of modes of the van Kampen type which seem not to arise within the framework of linear perturbation 
theory. Section 4 contains a very brief discussion of fully nonlinear oscillations of stellar systems, and the paper concludes in 
Section 5 with a summary and discussion of the important results of this work. 

Remarkably, many of the elements of what is fully worked out in the present investigation have been anticipated by 
Louis and Gerhard (1988), who asked the question 'Can galaxies oscillate?' and then constructed a self-consistent model of 
a spherically symmetric stellar system in a state of stationary oscillation. They constructed that model with the aid of a 
suitably generalized version of Schwarzschild's (1979) method for the solution of the fundamental integral equation describing 
the self-consistency of a stellar system. That paper provides an important example of a fully nonlinear (albeit numerical) 
model of a galaxy in a state of stationary oscillation. The work of Louis and Gerhard and its relationship to the present 
investigation are the main subjects of Section 4. It is relevant to the motivation for the present investigation to note that 
Louis and Gerhard (1988) begin their paper with a brief review of observational evidence that states of equilibrium may not 
be the only 'normal' states of galaxies. They suggest that some observed galaxies could be in states of stationary oscillation. 

The generality of the present investigation should be emphasized. In the first place, what follows establishes for finite, 
inhomogeneous stellar systems many results which have been well established for infinite, homogeneous plasmas and stellar 
systems. Moreover, these results apply to a wide class of stellar systems. The only significant restriction is that we consider 
the stationary oscillations of a system in which the motions of stars in the unperturbed gravitational potential are separable. 
For the sake of clarity and specificity, we include appendices in which the general analysis described in Sections 2 and 3 is 
worked out explicitly for plane waves in an infinite, homogeneous system of stars. 



2 MODES OF OSCILLATION OF THE VAN KAMPEN TYPE 

Infinitesimal perturbations in a galaxy are governed by the collisionless Boltzmann equation in the linearized form 



where the distribution function fo(x,v) represents the unperturbed density of stars in the six-dimensional phase space of 
a single star, Vo(x) denotes the unperturbed gravitational potential of the system, and fi(x,v,t) and Vi(x,t) denote the 
perturbations of the distribution function and the potential, respectively. Here, x and v denote the position and velocity of a 
star, respectively, and t denotes the time. The density of the perturbation 
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where m» is the mass of a star, is the source of the potential Vi(x,t). In other words, 

V 1 (x,t) = -G f pj^dx', (3) 
Jv \ x x \ 

where the integration extends over the volume V of the system. 

The normal modes of the system are described by solutions of equations (l)-(3) for fi(x, v,t), pi(x,t), and Vi(x,t) with 
a time dependence exp(— iat), where the characteristic frequency a is a constant. 

2.1 Solution for fi(x,v,t) in terms of action-angle variables 

We assume that the motion of a star in the unperturbed potential Vo(x) is separable. Accordingly, with the aid of a suitable 
canonical transformation (x, v) — > (w, J), we introduce a set of angle variables w = (wi, W2, W3) and conjugate action variables 
J = (Ji, J2, J3) which describe the motion of a star in the unperturbed system. When expressed as functions of the action-angle 
variables, functions of x and v are quasi-periodic functions of the angle variables w with periods 2n. 

The actions J are isolating integrals of the motion of a star in the gravitational field of the unperturbed system. Therefore, 
the unperturbed distribution function can be expressed as a function fo(x,v) = /o(J) of the actions in accordance with the 
theorem of Jeans. 

Transformed to action-angle variables, the linearized, collisionless Boltzmann equation now reads 
9/1 +u(J) dfl - 9Vl dfo (A) 

-m +u;{J) -d^-^-^j' (4) 

where the quantities w(J) = (loi,lo 2 ,lj 3 ) are the frequencies of the unperturbed motion of a star. In other words, w(J) = 
dH /dJ, where H (J) is the Hamiltonian which governs the unperturbed motion. 

In order to solve equation (4), we express the perturbations of the distribution function and the gravitational potential 
as trigonometric series in the angle variables of the forms 

fi(x, v, t) = fi(w, J,t) = /i( n > J' a ) cxp(in • w - iat) (5) 

n 

and 

Vi(x, t) = Vi(w, J, t) = ^i( n i J i a ) exp(in ■ w — iat), (6) 

n 

respectively, where the sums extend over sets of integers n — (m , 712, n.3). In this representation, trigonometric components 
of equation (4) separate, and the coefficients fi(n,J,a) and Vi(n, J, a) satisfy the equation 

-i[ff-n-«(J)]/i(n,J,«T) = in - ^Vi(n,J,a). (7) 

The solution of equation (7) for fi(n,J,a) must allow the integration of quantities involving fi(x,v,t) over the entire 
region of the phase space that is accessible to stars in the unperturbed system. In particular, /i(n, J, a) must be defined for 
this purpose at values of J such that a — n ■ = 0. For each set of integers n, that condition defines a surface of two 

dimensions in the three-dimensional space of the actions J. Let two variables (wi, V02) = to, say, label points on that surface. 
Then the values of the actions at a given point vj on the surface are given by the functions J = J r(ti, zu, a) (say). 

The solution for /i(n, J, a) is a generalized function of the form 

f 1 (n,J,a) = -P V ^ n > J >°l . n-ffe+ / A(n, w, a)8[J - J R (n, w, <r)]dw, (8) 

where P signifies that the integral of a quantity must be interpreted as the Cauchy principal value, A(n, vj, a) is an arbitrary 
function, and 8 denotes Dirac's delta function. The integral on the right-hand side of equation (8) extends over the entire surface 
in the action space on which a — n-u>( J) = 0. Equation (8) expresses /1 (n, J, a) as a sum of a particular solution of equation (7) 
and a general solution of the reduced equation [a — n • w(J)]/i(n, J, a) = 0. In particular, A(n, va, <r)5[J — J n(n, zu, a)]dzu 
is the contribution to the solution of the reduced equation (for fixed values of n and a) at values of the actions in the 
neighborhood of the values J = Jr(ti, zu, a). The integration over zu adds such contributions over all points J on the surface 
a — n ■ w(J) = 0. 

Equation (8) clearly satisfies equation (7) where er — n-w(J) 7^ and J / Jn(n,zu,a). Where a — n-w(J) = 0, equation 
(8) also satisfies equation (7) in the sense that 

{a - n ■ u[J R (n,zu,a)]}fi[n,J R (n,zu,a),a] 



lim 

J — *J^(n,w,(j) 



Uv-n-u(J)]\- Vl( > n > J >°l n .^Mp] + f \(n,zu',a)[a-n-u(J)]5[J-J R (n,zu',a)]dzu'} 
L I a — n ■ lj\J j OJ J J J 



(9) 



d j Vi(n,J,a) 
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Here the integral over vj vanishes in the limit, because 

a - n ■ w(J) — ► -[J - J R (n,T&',a)] ■ {-^[n ■ as J -> J R (n, vj' , a). (10) 

L OJ J J=J H (n,ro',<r) 

Substituting from equations (8) for the coefficients fi(n, J, a) in equation (5) we obtain the solution for the perturbation 
of the distribution function in the form 

fi(x,v,t) =fx(w,J,t) 

= — P ^("'^' ^l w . CX p(i n . _ -|_ I dzu\(n, zu, a)5[J — Jr(ti, tu, a)] exp(in ■ w — iat). 

It is in the solution of equation (7) that the present formulation of the matrix method differs from the formulation by 
Kalnajs (1977). Where the solution by Kalnajs is based on Landau's (1946) initial- value analysis of plasma oscillations, the 
present solution is based on van Kampen's (1955) normal-mode analysis. The difference is manifest in the interpretation of 
equations (8) and (11) in terms of Cauchy principal values of integrals and in the appearance in those equations of Dirac's 
delta function. Thus, for example, equation (11) is to be contrasted with equation (10) in Section 4 of the Appendix of Fridman 
& Polyachenko (1984) or with equations (9.5) and (9.6) in Palmer (1994). 



2.2 A matrix method for the solution of the characteristic value problem governing normal modes 

Substituting from equation (11) for fi(x,v,t) in equation (2) and rearranging terms in the resulting equation, we obtain 



i ,\ , D fj Vi(n,J,a) dfo ,. . ., 

Pi(x,t) + > m*P / dv r . n ■ —— exp(m • w — lot) 

^ / a — n ■ u>(J) dJ 

n 

= ^""^ m* J dv J dvj\(n, tu, a)8[J — Jr{i%, vj, a)] exp(in ■ w — iat). 



(12) 



Equation (12) is a condition of self-consistency to be satisfied by perturbations of the system, inasmuch as the density pi (a;, t) 
must be the source of the gravitational potential Vi(x,t). In other words, the characteristic value problem governing the 
normal modes of the system now reduces to equations (3) and (12). We can solve the characteristic value problem with the 
aid of a modified version of the matrix method of Kalnajs (1977). 

Let pi(x,/3) and Vi(x,f3) be members of a biorthonormal set of densities and potentials, in the sense of Clutton-Brock 
(1972) and Kalnajs (1976), where the constants /3 are sets of parameters which label the members of that basis set, and each 
density pi(x, (3) is the source of the corresponding potential Vi(x, (3). A biorthonormal set of densities and potentials satisfies 
the conditions 

f V 1 *(x,ct)p 1 (x,/3)dx = -5(a,l3), (13) 

where the integration extends over the volume V of the configuration, and 5(a,f3) is the Kronecker delta. The minus sign 
on the right-hand side of equation (13) is required by the physical interpretation of the integral on the left-hand side as 
the gravitational potential energy of the mass distribution described by the density pi(x,/3) in the gravitational potential 
Vi(x, a). When (3 = a, the integral represents the self-energy of the perturbation, which is negative. 

The matrix formulation of the characteristic value problem is based on representations of the perturbations of the density 
and the potential of the forms 

pi(x,t) = y^;u(/3, a)pi(x, /3) exp(-iat) and Vi(x,t) = p,{j3, a)Vi(x, j3) exp(-icrt), (14) 

respectively, where the coefficients p(/3, a) are constants to be determined. We shall need representations of the potentials 
Vi(x,f3) as trigonometric series in the angle variables of the form 

V 1 (x,/3) = V 1 (w,J,/3) = J2v 1 (n, J, (3) exp(in • w). (15) 

n 

It follows from the second of equations (14) that the coefficients in the trigonometric series on the right-hand sides of equations 
(6) and (15) must be related in the manner 

ft(n,J,ff) = 5^M/3,<7)Vi(n,J,/3). (16) 

D 

We obtain a matrix representation of the condition of self-consistency by multiplying each term in equation (12) by 
V{(x,a) and integrating the product over the volume of the configuration space accessible to stars in the unperturbed 
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system. In the reduction of the resulting equation, we transform integrals over the phase space from integrals over x and 
v to integrals over w and J. Inasmuch as this is a canonical transformation of coordinates and momenta, the Jacobian 
determinant of the transformation is unity. Making use of equations (13)- (16), and noting that integrals of terms of the form 
exp[i(n — n') • w] over the angle variables vanish unless n! = n, where n and n' are sets of integers, we reduce the resulting 
equation to 



= — ^^(27r) 3 m» J dwVi* [n, Jr(u, vj, a), a]X(n, z*7, a), 



(17) 



where we have suppressed a common factor exp(— iat). 



2.3 Normal modes 

Equation (17) represents a system of algebraic equations of the form 

5^M(ff,a,/3)Ai(/3,<r)=A(<7,a) (18) 



for the determination of the coefficients /i(/3,<r), where 

M(a, a, ffl = J(a, /3) - ^>^P J * ^ 0) „ . ^ dJ (19) 

and 

A(u, a) = — ^^(27r) 3 m* / Vi* [n, Jr(ti, vj, a), a]X(n, vj, a)dzj. (20) 

n ^ 

2. 5. J T/ie continuous spectrum of modes 

As they stand, equations (18)-(20) provide a matrix representation of the characteristic value problem governing modes with 
a continuous spectrum of real frequencies. For an assigned value of a in the continuous spectrum, equations (18) are an 
inhomogeneous system of linear equations governing the coefficients (j,(/3,cr). In general, those equations admit of a solution 
only if det|M(<r,a,/3)| + 0. 

Without loss of generality, we can normalize the perturbation in such a way that 



V 1 *{x,t) Pl (x,t)dx = -l, (21) 

v 

inasmuch as the gravitational potential energy of the perturbation must be negative. In virtue of equations (13) and (14), this 
normalization reduces to 

5>09,<7)| 2 = 1. (22) 



For the assigned value of a, equation (22) imposes one constraint on the quantities A(cr, a) and, hence, on the functions 
X(n,xu,a) (see equations [20]). Inasmuch as the functions X(n,zj,a) are arbitrary, apart from that constraint, there is no 
characteristic equation for the determination of the value of a. Moreover, the modes belonging to an assigned value of a in 
the continuous spectrum are highly degenerate. 

The physical meaning of the degeneracy of the modes in the continuous spectrum is the following. The functions X(n, vj, a) 
determine the perturbation fi(x,v,t) on the surfaces in the action space on which a — n ■ us(J) = and the unperturbed 
stellar orbits are resonant. The arbitrariness of the functions X(n,zu,a) represents a freedom to populate the resonant orbits 
in different ways in order to maintain the selfconsistency of the perturbation. 

One particular construction of modes of the van Kampen type can be achieved by choosing A(n, va, a) in the manner 

A(n, vj, a) = S(n, n )<5(tu - -n7 )A(n , vjq, a), (23) 

where no and tzto are fixed constants, 5(n, no) denotes a Kronecker delta, and S(zu — vuo) denotes Dirac's delta function. 
With this choice, we populate only the resonance n = no with stars, and, on the resonant surface a — no ■ w(J) = in the 
action space, we include only stars whose actions have the values J = Jr(uo, zjo, <t). In this case, equation (20) reduces to 

A(cr, a) = -(27r) 3 m» Vi* [no, Jjj(n , wo, o-),a]X(n , vjq, e). (24) 
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Therefore, all of the coefficients /i(/3, a) in the solution of equation (18) are proportional to A(no, vjq, a), ft follows that the 
normalization of the coefficients specified by equation (22) determines the value of |A(no, roo, c) | 2 - 

In Section 3, we shall find that nonlinear effects impose a constraint on the frequencies of the modes of the van Kampen 
type. Specifically, we shall find that the spectrum of frequencies is bounded by the frequency of a mode of the 'Vlasov type.' 

2.3.2 Modes of the Vlasov type 

At exceptional frequencies in the continuous spectrum, it can happen that det|M(er, a, (3) — 0. This is essentially a charac- 
teristic equation for the determination of such exceptional values of a. Under these conditions, a solution for the coefficients 
fi((3,a) exists if and only if A(a, a) = in equation (18). In other words the quantities X(n,zu,a) must all vanish. Never- 
theless, the integrals that appear in the definition of M(a,a,(3) in equation (19) must be interpreted as Cauchy principal 
values. In this case, the characteristic equation det|M(<r, a, f3)\ = is the counterpart of the dispersion relation proposed by 
Vlasov (1945) for plasma oscillations. This is a mode in the continuous spectrum in which stars on resonant orbits do not 
contribute to the self-consistency of the perturbation. The interpretation of Vlasov's dispersion relation in terms of an absence 
of particles trapped in a plasma wave has been discussed by Bohm and Gross (1949) and Jackson (1960). The interpretation 
of modes of the Vlasov type is illuminated by the discussion of their slightly nonlinear counterparts at the end of Section 3 
below. 

2.3.3 The discrete spectrum of modes 

In the cases of modes with complex frequencies and of certain modes with real frequencies, a — n ■ u>{J) ^ for all values of 
n and all values J accessible to stars in the unperturbed system. In such cases, there are no terms in the sums over integrals 
involving the functions X(n,zu,a) on the right-hand sides of equations (12), (17), and (20). In other words, the quantities 
A(cr, a) vanish identically, and equations (18) reduce to 



In this case, it is not necessary to interpret the integrals on the right-hand side of equation (19) as Cauchy principal values. 
Solutions of equations (25) for the coefficients n(f3,a) exist if and only if det|M(cr, a, f3)\ = 0. This condition provides a 
characteristic equation for the determination of the frequencies a of the discrete modes. 

For growing modes, the present version of the matrix method is equivalent to the version of Kalnajs, because the 
prescription implied by equation (11) for the integration of quantities involving the perturbation of the distribution function 
over the phase space coincides with Landau's prescription. Accordingly, the two versions of the matrix method give the same 
solutions for modes of instability. 

3 SLIGHTLY NONLINEAR OSCILLATIONS OF A GALAXY 

In the theory of modes of the van Kampen type described in the preceding section, we have attributed singularities in the 
perturbation of the distribution function to contributions of stars on resonant orbits. In the case of plasma oscillations, such an 
interpretation of the singularities in the van Kampen modes has been verified by Bernstein, et al. (1957) in their investigation 
of an exact model of nonlinear plasma waves (see also Stix 1992). 

In this section, we investigate stationary oscillations of galaxies in cases in which the amplitudes of the oscillations are 
small but finite. As in the investigation of Bernstein, et al. (1957), we find that such oscillations can be represented in terms 
of series in half-integral powers of the amplitude of the perturbation of the gravitational potential. We solve the governing 
equations through terms of the first order in that amplitude. The solution requires an explicit treatment of the contributions 
of resonant orbits to the density of stars in the six-dimensional phase space. This treatment of stationary oscillations confirms 
and clarifies the interpretation described above of the singularities in the modes of the van Kampen type. 

3.1 Formulation of the problem and reduction in terms of action-angle variables 

Consider an oscillating galaxy in which the density and the gravitational potential are of the forms 

p(x,t) = po{x) + ]-epi(x) exp(-krf) + \epl (x) cxp(iat) (26) 



respectively, where a is a constant, real frequency and the small, dimensionless parameter e is a positive, real number. We 
require that the oscillation be self-consistent in the sense that the density p(x, t) is the source of the potential V(x, t). Evidently, 




(25) 



and 

V(x,t) = Vo(x) + ^-eVi(x) cxp(-iat) + \eV?{x) exp(icrt), 



(27) 
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the density po(x) must be the source of the gravitational potential Vo(x). We shall find that po(x) and Vq(x) characterize 
an unperturbed state of the galaxy. Let fo(x,v) denote the distribution function of the galaxy in that unperturbed state. In 
equations (26) and (27), we represent p(x, t) and V(x, t), respectively, as real quantities, expressed, for the sake of convenience 
in what follows, in terms of a complex density pi(x) and a complex potential Vi(x), respectively. The requirement of self- 
consistency also implies that the density pi{x) must be the source of the potential Vi(x). 

The goal in what follows is to construct the distribution function f(x, v, t) of the galaxy in the state of oscillation described 
by equations (26) and (27). In order to accomplish this, we investigate the motions of stars in the gravitational potential 
described by equation (27), we construct a set of integrals of the motion, and, with the aid of the theorem of Jeans, we express 
f(x,v,t) as a function of those integrals. We carry out this program under conditions in which the oscillatory perturbation 
of the system is small but finite, and, in the construction of the distribution function, we evaluate all quantities involved 
consistently through terms of the first order in the parameter e. We have written the expressions for p(x, t) and V(x, t) as we 
have in terms of e in order to facilitate the ordering of terms in the series that arise in this process. 

As in Section 2, we consider a system in which the motion of a star in the unperturbed potential Vo(x) is separable, and 
we introduce the actions J and angles w of the unperturbed motion as canonical momenta and coordinates, respectively. The 
Hamiltonian is expressible as a function H = Hq{J) of the actions alone. It follows from the canonical equations of motion in 
that case that the actions J are integrals of the unperturbed motion of a star and the frequencies of the unperturbed motion 
are given by the relations u(J) = 8H /dJ. 

We now express the function V\ (x) as a function of the action-angle variables and particularly as a trigonometric series 
in the angle variables. Thus, we write 

Vi (as) = V x (w, J) = Vi (n, J) exp(in • w), (28) 

n 

where the sum extends over sets of integers n = (ni,n2,n.3). The motion of a star in the potential V(x,t) is governed by the 
Hamiltonian 

H(w,J,t) = Ho (J) + ^eVi(w,J) exp(-io-t) + ^eVi(w, J) exp(iat) 

1 , ( 29 ) 

= Hq(J) + y [Vi(n, J) exp(in ■ w — iat) + V 1 (n, J) cxp(— in ■ w + iat)]. 

n 



3.2 Construction of the integrals of the motion for a resonant orbit 

Consider now the motion of a star in the neighborhood of the resonance at which a — m ■ u>( J) = 0, where m = (mi,m2,m3) 
is a particular set of three integers. We can construct the integrals of the motion through order e in this case with the aid of 
two successive canonical transformations. This procedure is essentially von Zeipel's method for the construction of a canonical 
perturbation theory (von Zeipel 1916, Contopoulos 1975, Goldstein 1980). 

The first of these transformations replaces the actions J and angles w with new momenta i" = (ii, 12,1s) and new 
coordinates tp — (tj)i,'^2,'4>z) with the aid of the generating function 

" Vi(n, I) exp(in • w — iat) Vi(n, I) cxp(— in • w + iat)' 



F(I,w,t) = I w 



a — n ■ w(J) 



a — n ■ w(J) 



a function of the new momenta and the old coordinates. The transformation equations are 



_ dF T 1 

and 

8F i r 



"nVi(n, 7) exp(in ■ w — iat) nV 1 *(n, I) exp(— in ■ w + iat) 



a — n ■ w(J) 



a — n ■ w(J) 



(30) 



(31) 



Ud / V 1 (n,I) 
\ldl\a -n-w(I) 



cxp(in • w — iat) 



and the new Hamiltonian is 

K(ip, I, t) = H(w, J,t) + ^= H (I) + ieVi (m, I) cxp(im ■ V - iot) + ^eV*(m, I) exp(-im • i/> + iat) + 0(e 2 ) (33) 

(Goldstein 1980). In equation (33), we have eliminated the old angles w in favor of the new angles ip with the aid of equation 
(32). This canonical transformation removes the nonresonant terms of order e from the perturbed Hamiltonian. 

The second canonical transformation introduces final momenta G = (Gr,G2,Gi) and final coordinates 9 = {Or, 62,63), 
and it is designed to make the quantity m ■ ip — at & canonical coordinate. The generating function is 



F(G, ip, t) = G R (m ■ ij> - at) + G 2 i> 2 + G 3 V>3, 



(34) 



8 P. O. Vandervoort 



again a function of the new momenta and the old coordinates. The transformation equations are 

dF dF OF 

° R = ~^FT- =m-il> -at, 6 2 = -tttt = 1P2, and 63 = — —=^3 (35) 

oGr 0G2 0G3 

and 

dF dF dF 

Ii = -^r— = tuiGr, h — tt— =m,2GR + G2, and h — tt— — ttisGr + Gz, (36) 

ayji aip2 0^3 

and the final Hamiltonian is 

H r (9r, Gr, G ± ) = K(rP, /, t) + 7? 

1 l (37) 

= H (mG R + Gx) - oGr + -eVi(m, mG R + Gx) exp(iO R ) + -eV*(m, mG R + Gx) exp(-i0 H ) + 0(e 2 ) 

(Goldstein 1980), where we have defined Gx = (0, G2, G3) and we have abbreviated equations (36) to I = mGR. + G±. 

We observe that the new momenta G2 and G3 are integrals of the motion, inasmuch as the Hamiltonian Hr(0r, Gr, Gx) 
does not depend on the coordinates 82 and 63. Thus Hr(8r, Gr, Gx) is effectively a Hamiltonian governing a motion in one 
degree of freedom in the canonical variables (Or, Gr). The Hamiltonian Hr(9r, Gr, Gx) itself is a third integral of the motion 
inasmuch as it does not depend explicitly on the time. 



3.3 Reduction of the perturbed Hamiltonian to a pendulum model 

Let Go be the 'resonance value' of Gr satisfying the condition 

m • w(mG + Gx) - <r = (38) 

for assigned values of the components G2 and G3 of G±. We expect, and we shall find, that Gr — Go = 0(e 1//2 ) for stars on 
resonant orbits. In virtue of equation (38), an expansion of the Hamiltonian Hr(0r, Gr, G±) in powers of Gr — Go accordingly 
yields, in place of equation (37), 



Hr(0r, Gr, Gx) =H (mG + Gx) - aG + \{-^r\m- ■ w(mG + Gx)]}(G fl - Go) 



+ ieVi(m, mGo + Gx) cxp(i^) + \eV?(m, mG + G ± ) exp{-W R ) + Q(e 3/2 ). 



Let 



d 



[m ■ w(mGo + G±)] 



?? = signj^_-[m-w(mGo + Gx)]j and n= 

and write the function Vi(m, mGo + G±) in terms of its amplitude and phase in the manner 
Vi(m,mG + G ± ) = | Vi (m, mG + G±)\ exp(i0). 
Then, equation (39) reduces to 



H r (0r, Gr, Gx) = Ho (mG + Gx) - aG„ + r. 



-u(G R - G Q ) 2 + r)A Vi(m, mGo + G ± )\cos(0r + </>)+ Q(e i/2 ). 



(39) 

(40) 
(41) 
(42) 



This Hamiltonian has a stable fixed point at (Gr,0r) — (Go, So), where the value of 0o is such that cos(#o + 4>) = ~ r l- 
Inasmuch as the value of rj is either 1 or —1 (see equations [40]), it follows that an alternative form of equation (42) is 

Hr(0r, Gr, Gx) = H (mGo + G ± ) - aG„ + 77 [^(G fl - G f - e|Vi(m, mGo + G ± )\cos(0r - 6> )1 + 0(e 3/2 ). (43) 

This is essentially the Hamiltonian of a physical pendulum in the canonical momentum AG = Gr — Go and the conjugate 
coordinate AO — Or — On. The separatrix of the system, for an assigned value of G±, satisfies the equation 

\n{G R - Go) 2 - e|Vi(m,mGo + Gx)| cos(^ - 6> ) = t\Vi{m,mG + Gx)|, (44) 

where we are now ignoring terms of orders e 3 ^ 2 and higher in the Hamiltonian. In other words, the region of the phase space 
in which the angle Or librates about the value and the motion is resonant is denned by the inequality G_ < Gr < G + , 
where 



G ± = Go ± e 



1/2 



-|Vi(m,mG + Gx)| 



1/2 



[1 + COS(0K - 60)] 



1/2 



(45) 



This definition of the resonance region confirms the expectation described above that Gr, — Go = 0(e 1 ^ 2 ) for stars in the 
resonance region. In other words, the size of the resonance region is of order e 1 / 2 . Outside the resonance region, the angle Or 
precesses, and the motion is essentially nonresonant. 
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The resonance region described by equation (43) for the Hamiltonian and equation (44) for the separatrix consists of 
one or more islands in a closed chain in the phase space. The island or chain of islands circulates in the phase space. The 
guiding centre for the circulation of each island is a periodic orbit whose motion is given by J — tuGr + G± + O(e) and 
m ■ w — at + O(e) = 9 Q . 



3.4 A matrix representation of a self-consistent oscillation of a galaxy 

We are now ready to make use of the theorem of Jeans in order to construct the distribution function f(x, v, t) underlying the 
density distribution and the gravitational potential described by equations (26) and (27), respectively. We shall construct the 
distribution function in the non-resonance and resonance regions of the phase space separately. Moreover, we shall consider 
a system in which only a single resonance, defined by the set of integers m — (mi,m2,m3), is populated with stars. Finally, 
we shall develop the distribution function as a series in powers of e, working consistently to an order of approximation such 
that the fundamental relation 



m* J f(x, v, t)dv — p(x, t) (46) 

is satisfied through terms of order e. 

In this connection, we recall that we have identified the functions po(x) and Vo(x) in equations (26) and (27) as the density 
and gravitational potential, respectively, that characterize an unperturbed state of the galaxy. According to the theorem of 
Jeans, the distribution function fo{x, v) for the galaxy in that state can be expressed as fo{J), a function of the unperturbed 
action integrals. The self-consistency of the equilibrium of the galaxy in that state requires that 

m * J fo(x, v)dv = m* J fo(J)dv — po(x). (47) 
3.4-1 Non-resonant stars 

We first construct the distribution function for non-resonant stars, that is, for stars outside the region of the phase space 
defined by the inequality G_ < Gr < G+, where G± is defined by equation (45). For the construction of the integrals of 
the motion required for this purpose, we modify slightly the canonical transformation described by equations (30)-(32). In 
the trigonometric series in those equations, we now include the terms n — m which were previously omitted. In place of 
equation (33), the equation for the new Hamiltonian is now K(xp, I, t) — Hq(I) + 0(e 2 ). In this case, the canonical equations 
of motion imply that the new momenta I are integrals of the motion through terms of order e. According to the theorem of 
Jeans, the distribution function (for the non-resonant stars) can be expressed as a function of the integrals I. Inasmuch as 
the unperturbed distribution function of the system is of the form fo(x,v) — fo(J) described above, we write the perturbed 
distribution function for the non-resonant stars as 

, . . . . . 1 r^i( n > J) exp(in ■ w - iat) V* (n, J) exp(— in • w + iat) 1 df ~. 2 s /,„•> 

f N R(x,v,t) = fo(I) = fo(J) - -e > — — -f- + -f- n- — +0(e ) (48) 

2 ^-^ L a — n ■ u>(J) a — n ■ u>(J) J oJ 

n 

(see eq. [31], but recall that we now include the terms n = m in the trigonometric series there). 

The physical picture underlying this choice for the non-resonant distribution function is that the non-resonant orbits are 
orbits which are not altered qualitatively by the perturbation. In particular, the perturbed angle variables ip precess as do 
the unperturbed angle variables w. Thus, although the non-resonant orbits are perturbed, the dependence of the distribution 
function on the integrals of the motion is not altered for the non-resonant orbits. 

The contribution of the non-resonant stars to the density of the system is 



p NR (x,t) = m, / fo(I)dv 

J NR 

[ fo(J)dv-m*p[ L y^l V 1 (n,J)eMin-w-iat) + 

J NR. JnR 2 ^ L a - n ■ ^( J ) 



Vi (n, J) exp(in • w — iat) Vi (n, J) exp(— in ■ w + lat) 



■«(J) 



§^ + 0(^), 



where the integrals include only regions of the phase space in which the stellar orbits are non-resonant, and P signifies that 
integrals of functions containing resonance denominators are to be interpreted as Cauchy principal values. We reduce equation 
(49) to the more useful form 



p NR (x,t) = po(x) - m* / f {J)dv 



J) exp(in ■ w — iat) V{ (n, J) exp(— in ■ w + iat) 



n ■ u>(J) a — n ■ u>(J) 



n '§j dU + ° (e3/2) ' 



(50) 
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where we have rewritten the integral of fo(J) with the aid of equation (47) and we have made use of the label R in order to 
signify that an integral includes only the region of the phase space in which the orbits are resonant. In equation (50), we have 
also enlarged the region of integration of the last integral in equation (49) to include the region of the phase space containing 
resonant orbits. Inasmuch as the size of the resonant region is of order e 1 ^ 2 , the result is to change the integral by a quantity 
of order e 3 / 2 . That is beyond the order of approximation required here. 

The interpretation of integrals in equations (49) and (50) as Cauchy principal values now applies at the resonance n = m 
as well as at all other resonances. The reduction of equation (49) to equation (50) shows that, in general, the interpretation of 
the integrals as Cauchy principal values neglects the finite sizes of the resonance regions in the phase space and consequently 
neglects corrections of order e 3 ^ 2 to the density. 



3.4-2 Resonant stars 

We have found in Section 3.2 that the Hamiltonian Hr(9r,Gr,G±) and the components Gi and Gz of G± are integrals of 
the motion in the resonance region of the phase space. According to the theorem of Jeans, the distribution function for the 
resonant stars can be written as 

f R (x,v,t) = f R (H R ,G±), (51) 

where the function fR(HR,G±) remains to the determined (see Section 3.5). 
The contribution of the resonant stars to the density is 

p R (x,t) = m* / f R (H R ,G±)dv. (52) 

J R 



3.4-3 A matrix representation of self-consistency 

The condition that the oscillation of the galaxy be self-consistent requires that the sum of the densities pnr(x, t) and Pr(x, t) 
be equal to the 'imposed' density p(x,t). Combining equations (26), (50), and (52), canceling the term po(x) on the two sides 
of the resulting equation, and rearranging terms in what is left of the resulting equation, we obtain 



i* J fR(H R ,G±)dv - m» J f (J)dv = ^epi(x) exp(-iat) + ^ep{(x) exp(iat) 

f 1 T^i( n ! J) exp(in • w — iat) V*(n, J) exp(— in • w + iat) 



n ■ u>(J) a — n ■ w(J) 



§^d. + 0( e 3 / 2 ). 



(53) 



Equation (53), which expresses the condition of self-consistency for a slightly nonlinear oscillation of the system, is to be 
compared (apart from terms of orders e 3 ^ 2 and higher) with equation (12), which expresses the condition of self-consistency 
for an infinitesimal perturbation of the system. In particular, the right-hand side of equation (53) is essentially the real part of 
the left-hand side of equation (12). Thus, the essential difference between the two equations lies in the different representations 
of the contributions of the resonant stars on the right-hand side of equation (12) and on the left-hand side of equation (53). 

Beginning with equation (53), we can now formulate the matrix method for the determination of the functions p±(x), 
Vi(x), and f R (H R , G±). In this formulation, we represent pi(x) and Vi(x) in terms of the biorthonormal set of densities and 
potentials introduced in Section 2.2. In particular, we now let equations (13)-(16) apply to pi(x,t) — pi(x) exp(— iat) and 
Vi(x,t) = Vi(x) exp(— iat). 

We next multiply each term in equation (53) by Vi(x, a) exp(kri), integrate the result over the configuration space, and 
then average the result over one period 27t/<t of the oscillation. The reduction here of the terms appearing on the right-hand 
side of equation (53) follows closely the similar reduction in Section 2.2 of terms appearing on the left-hand side of equation 
(12) and leading to the terms on the left-hand side of equation (17). In other words, the immediate result of these operations 
on the terms in equation (53) is 

/2n I <? {• 
dt exp(krf)m„ / dxdvVi {x, a)[f R (HR, G±) - fo( J)] 

(54) 



V 1 *(n,J,a)V 1 (n,J,f3) n dfo 
a — n ■ u>{J) dJ 



p(P,a) + 0(e i /' 2 ) 



On the left-hand side of equation (54), we transform the integral over the resonance region of the phase space to an 
integral over the action-angle variables J and w and we recall that the Jacobian determinant of the transformation is unity. 
Writing Vi(x,a) as a trigonometric series in the angle variables in accordance with equation (15), we obtain 



m„ / dxdvV{(x, a)[f R {H R , G±) - fo(J)] = ^ m« / dwdJV{(n, J, a) exp(-in ■ w)[f R (H R , G±) - fo(J)]- 
Jr Jr 



(55) 
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We next transform the integral over the actions J on the right-hand side of equation (55) to an integral over the momenta 
Gr, G2, and G3. For this purpose, we make use of the transformation equations (31), (32), (35), and (36), and we observe 
accordingly that 

J = mG R + G± + O(e) = mG + G± + m(G R - Go) + 0(e), and dj = mi[l + 0(e)]dGijdG 2 dG 3 . 

In the transformed integral, the integration over Gr is restricted to the resonance region G_ < Gr < G+, where G± is given 
by equation (45). Inasmuch as Ji = miCn + 0(e), the limits of integration over Gr depend on the algebraic sign of mi. Thus, 
if mi > 0, then the upper and lower limits of integration are G+ and G_, respectively. However, if mi < 0, then the upper 
and lower limits of integration are G_ and G+, respectively. In order to include both cases in equation (56) below, we replace 
midGfl with \mi\dGR and consistently adopt G+ and G_ as the upper and lower limits of integration, respectively. 
With the further observation that 



fo{J) =/o(mG + Gi) + 
and 



J-/o(mGo + Gi) (G B - Gn) -O(c) 



V*(n,J,a) = Vi* (n, mGo + G± , a) + 
we finally reduce equation (55) to 

, f dxdvV 1 *(x,a)[f R (HR,G±) - fo(J)] 
Jr. 



d 

-^7- Vi (n, mGo +Gi,a) 



(G fl - Go) + O(e) 



: ^^ m * J dt« cxp(— in • to) ^ 



dtyexp(-in-«;) / dG 2 dGsW(n,mGo + G ± , a)|mi| / [f R (H R , G±) - /o(mG + G±)]dG H [l + O(e)]. 



(56) 



Contributions of order e 1//2 vanish on the right-hand side of equation (56), because they involve integrals over Gr of odd 
functions of (Gr — Go) with limits of integration which are symmetric with respect to the value Gr = Go- 
In the next subsection, we shall seek a solution for fR(H R , G±) such that 

f°+ 1 1 

\mi\ / [f R (H R , Gi) - /o(mG + G ± )]dG R = -eA(ro,G±, Go) exp(i0K) + -eA* (m,G±, Go) exp(-i0«), (57) 
j g _ 

where A(m, G±, Go) is a function to be determined. In the remainder of this subsection, we shall complete the formulation of 
the matrix method. Accordingly, we now combine equations (56) and (57) in order to rewrite the left-hand side of equation 
(54). Inasmuch as exp(i^ij) = exp(im ■ w — \at)[l + 0(e)] (see eqs. [32] and [35]), we bring equation (54) to the form 



/< 



i(27r) 3 em. / dG .dG- \ V (m. m('V +G..tt)A(m.G i .G„) 



-e^[5(a,/3)-^( 27 r) m.P / dj ^T^J) " ' J M/3 ' a) + ° (e } ' 

j3 n 



(58) 



Equation (58), with the terms of orders e 3 ' 2 and higher suppressed, is an inhomogeneous, linear equation for the de- 
termination of the coefficients /x(/3, a). To order e, equation (58) is essentially the matrix representation of self-consistency 
obtained previously in equation (17) for modes of the van Kampen type. The most substantial difference between equations 
(17) and (58) is that the inhomogeneous terms on the right-hand side of equation (17) represent a sum of contributions of 
resonant stars, whereas the inhomogeneous term on the left-hand side of equation (58) represents the contribution of a single 
resonance. The remaining differences between the two equations are entirely superficial. The components G2 and G3 of G±, 
which label points on the surface a — m ■ u(mGo + G±) = in the action space, are equivalent to the components w\ and 
VJ2 of -UJ, which label points on the surfaces a — n ■ w( J) = 0. Likewise, the quantity mGo + Gi is equivalent to the function 
Jf?(n, vj, a). The dependence of A(m, G±, Go) on Go is implicitly a dependence on m, a, and Gi in virtue of equation (38). 

Equation (58) is readily generalized in order to include contributions of a number of resonances on the left-hand side. 
For this purpose, the analysis in Sections 3.2, 3.3 and 3.4.2 would be performed within each of the resonance regions that 
contain stars, and the results incorporated in equations (53)-(58). In each of the resonance regions included, the solution for 
the perturbation of the distribution function would be determined as is described in Section 3.5. 

The solution of equation (58) for the coefficients fj,(/3,a) can be obtained and the constraint on A(m, G±,Go) imposed 
along the lines described in Sections 2.2 and 2.3 in the case of equation (17). However, in the case of slightly nonlinear oscil- 
lations, we shall find in the next subsection that the solution of equation (57) for /r(Hr, Gi) imposes additional constraints 
on A(m, G±, Go) which do not arise in the linear theory of modes of the van Kampen type. 



12 P. O. Vandervoort 



3.5 Construction of the distribution function for the resonant stars 

We shall now solve equation (57) for f R (H R , G±). Let 

f R (H R ,G±) = fo(mGo + G±) + F R {E R ;m,G±,G ), (59) 
where F R (E R ; m, G±, Go) is the component of f R (H R , G±) which remains to be determined, we have defined an 'energy' 

E R = ri{H R (e R , G R , Gi) - [H {mGo + G±)- oG Q + »je|Vi(m, mG Q + G±)\]} 

1 o (60) 

= -M(AG) 2 - e|Vi(m,mG + G±)|[l + cos(A0)] 

(see eq.[43]), and we have introduced new canonical variables AG = G R — Go and AO = 6 R — 9o- Then we can rewrite equation 
(57) in the form 

fG + -G 

\mi\ / Fr(E r ; m, G±, Go)d(AG) = eRe[A(m, G±, G ) exp(i0 o )] cos(A0) - eIm[A(m, G±, G ) exp(i0 o )] sin(A0). (61) 

It follows from equations (45) and (60) that the integral on the left-hand side of equation (61) must be a function of cos(A#) 
and, hence, an even function of A9. Therefore, equation (61) can be satisfied only if 

Im[A(m, G±, Go) exp(i0 o )] = 0. (62) 
In other words, the quantity A(m, G±, Go) exp(i#o) must be real, and equation (61) reduces to 

rG + —G 

\mi\ F R {E R ;m,G ± ,Go)d{AG) = eX{m,G±,G o )exp(i0 o )cos{Ae). (63) 

Equation (63) can be transformed into an integral equation of the Abelian form. Making use of equation (60) in the inversion 
of equation (63), we obtain the solution 

»"* G - °°> = -(f) w? w °X^S> w <- " G ° + «m(-«-> - 1/1 - '«» 

Upon combining this result with equation (59), we obtain the distribution function for the resonant stars in the form 

MH R ,G ± ) = MmG + G ± )-(^) 1/2 [e| V, (rrz, rrzg + g^) | (-^) ~ _ 2 (-^) V^j , (65) 

\2J \m\ [K\V\(rn, mGo + G-j_)| 

The term /o(mGo + G±) on the right-hand side of equation (65) is of order unity, and it describes a uniform density of 
stars in the resonance region of the phase space. Moreover, the contribution of /o(mGo + G±) to the density of stars in the 
resonance region and the contribution of fo(J) to the density in the non-resonance region (see eq. [48]) combine in equation 
(46) to give the unperturbed density po{x) on the right-hand side of equation (26). 

The second term on the right-hand side of equation (65) is formally of order e 1,/2 , inasmuch as E R is of order e. That 
contribution to the distribution function in the resonance region and the contribution of the terms of order e on the right-hand 
side of equation (48) to the distribution function in the non-resonance region combine in equation (46) to give the sinusoidal 
contributions to the density on the right-hand side of equation (26). 

The second term on the right-hand side of equation (65) is mildly singular as E R — * 0. Therefore, the requirement that 
the total density of stars in the six-dimensional phase space must be non-negative implies that 

A(m, G ± , Go) exp(i0 o ) < 0, (66) 

where it will be remembered that A(m, G±, Go) exp(i#o) is a real quantity (see eq. [62]). 

The singularity of f R (H R , G±) at E R = implies that the resonant stars tend to concentrate near the separatrix defined 
by equation (44) (cf. eq.[60]). 

If A(m, Gx, Go) =0, then the oscillation described here is a slightly nonlinear counterpart of a mode of the Vlasov type. 
According to equation (58), resonant stars are not required to support the oscillation of the system. Moreover, the mode 
satisfies the characteristic equation det|M(a, a,0)\ = as described in Section 2.3.2. According to equation (65), the density 
of stars in the resonance region of the phase space is uniform in that case. 



3.6 A constraint on modes of the van Kampen type 



Inequality (66) constrains the possible frequencies of a stationary oscillation described by equation (58). In order to derive the 
constraint, we multiply each term in equation (58) by n*(a,a), sum over the values of a, and reduce the resulting equation 
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with the aid of equations (16), (22), and ((41). Finally, we make use of the condition cos(#o + (j>) — governing the value of 
9r — 6q at a stable fixed point of the Hamiltonian represented by equation (42) (see also equations [40]). We obtain 



It follows from inequality (66) that the product of rj and the left hand side of equation (67) must be negative or zero. Therefore, 



The equality holds in condition (68) in the case that X(m, G±,Go) = 0. In other words, condition (68) implies that the 
continuous spectrum of frequencies of the modes of the van Kampen type is bounded by the frequency of a mode of the Vlasov 
type. The continuous spectrum might be bounded either from above or from below, depending on the sign of r\ (see equations 
[40]) and on the structure of the integrand in condition (68). 

The constraint imposed by condition (68) on the frequencies of plane waves in an infinite homogeneous stellar system is 
described in Appendices B and C. 

4 NONLINEAR OSCILLATIONS OF A GALAXY 

In principle, the series in powers of e constructed in Section 3 for the representation of slightly nonlinear oscillations could be 
extended to include terms of higher orders in e. However, there does not appear to be any outstanding issue of principle that 
would require the investigation of terms of higher orders. Moreover, the construction of terms of higher orders in the series 
representation of slightly nonlinear oscillations of a galaxy would rapidly become complicated. Therefore, it seems preferable 
to find alternative ways to satisfy equations (26), (27), and (46) for fully nonlinear oscillations without resorting to the series 
solutions started in Section 3. 

For a particular model of a radial oscillation of a spherically symmetric galaxy, Louis and Gerhard (1988) have solved 
equations (26), (27), and (46) with the aid of a generalized version of Schwarzschild's (1979) numerical method for the 
construction of models of galaxies. Specifically, they have made use of Lucy's (1974) method in order to solve the discrete set 
of equations with which Schwarzschild replaces the integral equation (46) . 

In their discussion and interpretation of their results, Louis and Gerhard touch upon certain features of their model which 
are reproduced more generally in the present work. By investigating the arrangement of resonant orbits in the phase space 
and exhibiting the solution of Schwarzschild's equations for the distribution of stars on the resonant orbits, they demonstrate 
that the self-consistency of the oscillation is sustained by an appropriate population of resonant orbits with stars. In that 
connection, they present and analyse the pendulum model described in Section 3.3 above in the form appropriate to conditions 
of spherical symmetry. Their surfaces of section exhibit the arrangement of the resonance regions of the phase space as closed 
chains of islands as described at the end of Section 3.3 above. Finally, Louis and Gerhard remark that resonant 'orbits near 
the unstable closed resonant orbits support the perturbation.' This accords with the more general result found at the end of 
Section 3. 5 above that the resonant stars tend to concentrate near the separatrix defined by equation (44). 

5 DISCUSSION AND CONCLUDING REMARKS 

The work described in this paper generalizes the theory of van Kampen modes in plasmas to the gravitational case of 
modes of oscillation in inhomogeneous stellar systems, and it does so with the aid of a modified formulation of the matrix 
method of Kalnajs. The present construction of normal modes, particularly modes with continuous spectra of real frequencies, 
complements more conventional, initial- value treatments of small perturbations in galaxies. 

We have found that the connection between van Kampen modes and nonlinear plasma waves carries over to a connection in 
stellar dynamics between linear modes of stationary oscillation and slightly nonlinear, stationary oscillations. That connection 
also bridges the gap between the linear theory of small perturbations in stellar systems and nonlinear models of oscillating 
systems of the kind considered by Louis and Gerhard (1988). The slightly nonlinear theory of stationary oscillations described 
in Section 3 should provide an analytical foundation for more general and more extensive numerical investigations along the 
lines described by Louis and Gerhard. 

The treatment of stationary oscillations described in this paper is quite general, and it has wide applications. The principal 
restriction underlying the construction of linear oscillations in Section 2 and slightly nonlinear oscillations in Section 3 is that 
the stellar motions in the unperturbed gravitational field of the system must be separable. The KAM theorems imply that the 
series in powers of the small parameter e developed in Section 3 are a valid approximation, for sufficiently small values of e, 
in the sense that the perturbed motions of stars remain separable for most initial conditions. Thus, the theories of stationary 
oscillations of stellar systems described in this paper apply to spherically symmetric systems and to axisymmetric and triaxial 




n 




(68) 



n 
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Stackel systems. The KAM theorems also imply that the present representation of stationary oscillations in galaxies would 
apply, at least approximately, to systems in which the unperturbed stellar motions are only approximately separable. This 
would include applications to systems such as slowly-rotating, triaxial systems. 

On the other hand, for sufficiently large values of e, the overlap of resonances will lead to chaotic behavior of the relevant 
stellar orbits (Chirikov 1979). Such chaotic behavior would be expected to inhibit stationary oscillations at sufficiently large 
amplitudes. Inasmuch as orbits near the separatrices that form the boundaries of resonance regions would be the first to 
become chaotic, the onset of chaos in the oscillating systems considered here would probably be enhanced by the tendency of 
resonant stars to concentrate near the separatrices. 

Although the oscillations described in this paper are sinusoidal, it appears that more general periodic oscillations would 
arise in a more fully nonlinear treatment. Louis and Gerhard (1988) showed this to be the case for their model, but they argued 
that the relevant nonlinearities were small enough that they could approximate the oscillation of the model as sinusoidal. 

It is an important result of the present investigation that nonlinear considerations impose the constraint represented by 
inequality (68) on the continuous spectrum of modes of the van Kampen type and that a mode of the Vlasov type plays 
a special role in this connection. The mode of the Vlasov type is the limiting case of modes of the van Kampen type in 
which resonant stars do not contribute to the self-consistency of the oscillation. Moreover, according to inequality (68), the 
continuous spectrum of frequencies of modes of the van Kampen type is bounded by the frequency of the mode of the Vlasov 
type (see particularly the example described in the appendices of this paper). These results, which seem not to be recognized 
in the literature of plasma physics, raise an important issue regarding the connection between Landau's and van Kampen's 
treatments of plasma oscillations. The point of inequality (68) is that the excluded modes of the van Kampen type are 
unphysical in the sense that involve negative densities of resonant stars in the phase space. One might ask if unphysical van 
Kampen modes can be included in superpositions of modes that would represent solutions of the initial value problem of 
the Landau type. There is certainly no mathematical objection to the use of unphysical modes in a superpositon, and there 
can be no physical objection if the superposition preserves positive densities of stars everywhere in the phase space. Thus, 
inequality (68) serves mainly to identify those modes of the van Kampen type which, in isolation, may be regarded as physical 
perturbations. 

Stix (1992, see particularly Chapters 7 and 8) has emphasized that the connection between van Kampen modes and BGK 
waves validates the consideration of the van Kampen modes individually as distinct modes of oscillation in a plasma. The 
results described in this paper serve likewise to validate the consideration of stationary oscillations of galaxies. It remains, 
however, to identify mechanisms that could excite stationary oscillations of galaxies in nature. That is a question for future 
research. Nevertheless, one might speculate that, in the formation or growth of a galaxy through mergers, a stationary 
oscillation could be excited as an accreted fragment becomes trapped in a resonance. And one might ask if dynamical friction 
acting on an accreted fragment could contribute to such trapping. 
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APPENDIX A: VAN KAMPEN MODES IN AN INFINITE, HOMOGENEOUS SYSTEM 

Perturbations of a system in which the density is uniform and the distribution function fo(v) depends only on the stellar 
velocities provide an illustration of the matrix method described in the body of this paper. The model is the gravitational 
counterpart of the theory of plane waves in a homogeneous plasma (van Kampen 1955; Stix 1992). The evolution of small 
perturbations in an infinite, homogeneous stellar system has been treated as an initial-value problem by Lynden-Bell (1962) 
and Sweet (1963), and the complementary normal-mode analysis is discussed briefly by Binney and Tremainc (1987). 

In order to construct a suitable set of action-angle variables and a suitable biorthonormal set of densities and potentials, 
we employ a 'box normalization' and concentrate on the unit cell < Xi < L (i = 1, 2, 3). The action-angle variables of the 
unperturbed motion of a star are 

J — — and w = — - — . (Al) 

27T L 

The Hamiltonian and the frequencies of the unperturbed motion are 

Ho(v)= 1 -\vf= 1 -(^y\jf = H (J) and «,( J) = ^ = (£) *J = (§) v, (A2) 
respectively. The members of the basis set are of the form 

Ikl 2 

Vi(x, k) = Vi(k) cxp(ik ■ x) and pi{x, k) — pi (ft) cxp(ik ■ x) = — — q^( k ) exp(ik • x), (A3) 

where n — (ki, K2, K3) is a set of three integers and k = 2-kk/L. Thus the basis functions satisfy periodic boundary conditions 
on the unit cell. In the orthogonality relations, equations (13), the unit cell is the volume of integration. In particular, the 
normalization of these bases gives 

|^|Vi(K)| 2 L 3 = l. (A4) 

We must reduce equations (11) and (17) with the aid of equations (A1)-(A4). Upon writing equation (15) explicitly for 
this purpose, we find that 

V 1 (n,J,f3)=5(n,/3)V 1 (f3), (A5) 

where 5(n,/3) is a Kronecker delta. When we simplify equation (17) with the aid of equation (A5)), we find that the modes 
must be plane waves exp(ik ■ x — iat) — exp(i« ■ w — iat). Thus, in equations (11) and (17), we have fj,(a,a) = and 
A(a, o,ir) = 0(a/ #e) and, in virtue of equation (22), fi(n, a) = 1. This reduces each of the trigonometric series in equation 
(11) to a single term in exp(i/t ■ w). 

The condition a — n ■ u(J) = in the case that n = k reduces to a — k ■ v = 0. By resolving v into a component Vk 
parallel to k and a component v± perpendicular to k, we can reduce the conditions J — Jn(n,zu,a) to Vk = cr/|k| and 
v± — VJ. 

Equation (11) reduces to 

fi(x,v,t) — — P YUlII k . ^2. exp(ik ■ x — iai)+(^\ \(k,v±, a)8(vk — c/|k|) exp(ik • x — iat). (A6) 
cr — K - V OV \ L J 

Similarly, equation (17) reduces to 

1 ~^W ZP j dV a -k v k ' H = -(27r) 3 m. J dv ± X(K,v ± ,a). (A7) 

Equations (A6) and (A7) describe the gravitational analogue of van Kampen's solution for the electrostatic modes of oscillation 
of a homogeneous plasma (van Kampen 1955). 

For a mode in which the frequency a is real, equation (A7) fixes the value of the integral of X(k, v±,a) over v± on the 
right-hand side, and there is no dispersion relation. Of course, the exceptional case is one in which X(k, v±,a) vanishes, and 
equation (A7) reduces to the gravitational counterpart of Vlasov 's (1945) dispersion relation. 
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When |fe| < kj, where kj is the Jeans wave number, and a — k • v 7^ 0, equation (A7) reduces to 

47rGm* f 1 dfo „ , . 

1 -^k^7 du ^k^ k i7= ' (A8) 

a dispersion relation which admits unstable Jeans modes (Binney and Tremaine 1987). If the velocity distribution fo(v) is 
Maxwellian, then kj 2 — 12-ixGp/ {v 2 ) , where (v 2 ) is the mean-square velocity of a star in three dimensions. 



APPENDIX B: BGK WAVES IN AN INFINITE, HOMOGENEOUS SYSTEM 

In this appendix, we construct a gravitational counterpart of a BGK wave (Bernstein, et al. 1957) in an infinite, homogeneous 
stellar system as an example of a slightly nonlinear oscillation of the kind described in Section 3. 

As in the case of a van Kampen mode of the kind described in Appendix A, the perturbation of the system consists of a 
single plane wave. Symbols not defined in what follows have the meanings assigned in the body of the paper or in Appendix 
A. Thus the density of the system is of the form 

p(x, t) = po + ^epi(n) exp(ik ■ x — iat) + iep*(/«) cxp(— ik ■ x + iat), (Bl) 

where po is the constant density of the system in the absence of the wave, and the gravitational force per unit mass acting on 
a body is derived from the potential 

V(x,i) = ieVi(K;)exp(ik ■ x - iat) + ^eV*(n) exp(-ik ■ x + iat). (B2) 

Here, the wave is described in terms of a single member of the biorthonormal basis set of potentials and densities represented 
in equations (A3). The Hamiltonian is now of the form 

H = TjM 2 + ^Vi(k) cxp(ik ■ x - iat) + ^eV* (k) cxp(-ik ■ x + iat). (B3) 

In virtue of equations (Al) and (A2) the resonant condition reduces to a — k ■ v = in this case. The Hamiltonian 
contains no non-resonant terms. Therefore, the canonical perturbation theory described in Section 3.2 simplifies so that the 
transformation equations (Al), (31), (32), (35), and (36) reduce to 

vL 

— — J — I — kGr + G±, 

fB4) 

. . 2TTX2 , „ , 271-23 

Or — K ■ ip — at — K ■ w — at — k ■ x — at, O2 = ip2 = W2 = — ; — , and 6/3 = y>3 = W3 = — - — , 

and we are left with the Hamiltonian 

1 / 077- \ 2 1 1 

H r (8r, Gr, G±) = - (-) \kGr + Gx| 2 - aG R + -eVi(is) exp(i^) + -el? (is) exp(-i^) (B5) 

(see eqs. [37] and [A5] and note that we are now replacing m with k). 

In the reduction of the right-hand side of equation (B5) to the Hamiltonian of a pendulum model, we find that equations 
(40) reduce to r) = +1 and p, — \k\ 2 in virtue of the second of equations (A2) and the first of equations (B4). As in Appendix 
A, we introduce components Vk and v± of v parallel and perpendicular, respectively, to the wave vector k. With the aid of 
the first of equations (B4), we find that 



v k = kG R + ^-k ■ G± and v± = ^ 
kL L 



fe(fe-Gx) 

LrX 



(B6) 



fc 2 

Thus for the Hamiltonian of the pendulum model in the case at hand, we replace equation (43) with 

Hr(0r,Gr, Gx) = ^v 2 + 7}\ v ±\ 2 - & y + a ^~~j~tr~ + \( Vk _ u °) 2 _ e|Vi(/t)|cos[fc(a;fl - x )], (B7) 

where vo = a/k is the phase velocity of the wave. Here, we have also defined xr and xo by writing Or = kxR — at and 
#0 = kxo — at, respectively, where it will be remembered that 60 is the value of the angle variable 9r at the stable fixed point 
of the Hamiltonian. 

The distribution function for the stars trapped in the wave, i. e. , for the resonant stars, is now 
fR(H R ,G ± ) = M vo) - Q) V2 NM*)l(-^r 1/2 ~ 2(-E R ?'*\, (B8) 

where vo = (kvo/k) + v± and 

E R = i(w fe - v ) 2 - e|Vi(/t)|{l + cos[k(x R - x )]} (B9) 
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(see eqs.[60] and [65]). 

The distribution function for the untrapped stars, i. e. , for the nonresonant stars, is 



f N R(x,v,t) = f (v) - -e 



Vi(k) exp(ifc • x — iat) V*(n) exp(— ife ■ x + iat) 



+ 



k ■ 



dfo 
dv 



(BIO) 



a — k ■ x a — k ■ x 

(see eq.[48]). 

Finally, inasmuch as the perturbation of the potential consists of a single basis function Vi(#t) exp(ik • x — iat), equation 
(58) expressing the self-consistency of the wave reduces to 



(2-Kfem, / dG2dG 3 V 1 *(K)\(K,G±,G ) = 



47rGm* 



f^^k- 

J a — fe • v 



k 



dfo 
dv 



in virtue of equations (Al) and (A4). Equation (Bll) is to be compared with equation (A8). 
For a plane wave described by equation (Bll), inequality (68) reads. 



1 - 



4nGm 



IP [ dv 1 k .^l< . 

J a — k ■ v ov 



(Bll) 



(B12) 



Where the equality holds, this condition gives the gravitational counterpart of Vlasov's dispersion relation for a plasma wave. 
The manner in which the mode of the Vlasov type bounds the allowed, continuous spectrum of modes of the van Kampen 
type is illustrated in the appendix that follows. 



APPENDIX C: AN EXAMPLE OF AN ALLOWED SPECTRUM OF VAN KAMPEN MODES 

We define a one-dimensional velocity distribution 

FoM = / f (v )dv± (CI) 



for the unperturbed system, and we let 

F (v k ) = ^-(l-j 2 v 2 k ) (vl<j- 2 ) and F (v k ) = (v 2 k > j' 2 ), (C2) 
where j is a constant. 

For this model, the integral in inequality (B12) is readily evaluated and the inequality brought to the form 



faGA,-*L*- ta (i£l)]' (C3) 
where y = k/(aj) and, without loss of generality, we let a > 0. 

For an assigned value of y, the right-hand side of inequality (C3) sets an upper bound a max on the magnitude of a. Thus, 
o — <7 max at a wave number k — ajy. We have a = a max when the equality holds in conditions (B12) and (C3). In other 
words, Umax satisfies the gravitational counterpart of the Vlasov's (1945) dispersion relation. 

Figure Al is a plot of (J max against k for a plane wave in the model that we are considering. Here, we measure frequencies 
in the unit (QnGpo) 1 ^ 2 and wave numbers in the unit (6nGpoj 2 ) ll/2 ■ The frequencies in the continuous spectrum of the van 
Kampen modes which are allowed by inequality (C3) fill the region under the curve in Figure Al. Note that a max — + as 
k — > kj = (12-KGpoj 2 ) 1 ^ 2 , where kj (= 2 1 ^ 2 in the adopted system of units) is the Jeans wave number for the model. Evidently, 
the allowed modes of the van Kampen type occur in this model only at wave numbers such that k < kj and the system is 
gravitationally unstable. 



This paper has been produced using the Royal Astronomical Society/Blackwell Science TpX macros. 
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Figure Al. The allowed real frequencies in the continuous spectrum of modes of the van Kampen type in an infinite, homogeneous 
stellar system. For a plane wave, the curve represents the solution of a dispersion relation of the Vlasov (1945) type for the dependence of 
the frequency <r, measured in the unit (6-irGpo) 1 / 2 , on the wave number k, measured in the unit (Q-wGpoj 2 ) 1 / 2 . The allowed frequencies 
of the van Kampen modes fill the region under the curve representing the Vlasov mode. 
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